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Many eukaryotic cells chemotax, sensing and following chemical gradients. However, experiments 
have shown that even under conditions when single cells cannot chemotax, small clusters may still 
follow a gradient. This behavior has been observed in neural crest cells, in lymphocytes, and dur¬ 
ing border cell migration in Drosophila, but its origin remains puzzling. Here, we propose a new 
mechanism underlying this “collective guidance”, and study a model based on this mechanism both 
analytically and computationally. Our approach posits that the contact inhibition of locomotion 
(CIL), where cells polarize away from cell-cell contact, is regulated by the chemoattractant. Indi¬ 
vidual cells must measure the mean attractant value, but need not measure its gradient, to give rise 
to directional motility for a cell cluster. We present analytic formulas for how cluster velocity and 
chemotactic index depend on the number and organization of cells in the cluster. The presence of 
strong orientation effects provides a simple test for our theory of collective guidance. 


Cells often perform chemotaxis, detecting and mov¬ 
ing toward increasing concentrations of a chemoattrac¬ 
tant, to find nutrients or reach a targeted location. This 
is a fundamental aspect of biological processes from im¬ 
mune response to development. Many single eukaryotic 
cells sense gradients by measuring how a chemoattrac¬ 
tant varies over their length [T]; this is distinct from 
bacteria that measure chemoattractant over time [2]. In 
both, single cells are capable of net motion toward higher 
chemoattract ant. 

Recent measurements of how neural crest cells respond 
to the chemoattractant Sdfl suggest that single neural 
crest cells cannot chemotax effectively, but small clusters 
can [5]. A more recent report shows that at low gra¬ 
dients, clusters of lymphocytes also chemotax without 
corresponding single cell directional behavior; at higher 
gradients clusters actually move in the opposite direction 
to single cells [1]. In addition, late border cell migration 
in the Drosophila egg chamber may occur by a similar 
mechanism [Mg. These experiments strongly suggest 
that gradient sensing in a cluster of cells may be an emer¬ 
gent property of cell-cell interactions, rather than arising 
from amplifying a single cell’s biased motion; interest¬ 
ingly, some fish schools also display emergent gradient 
sensing |g. In fact, these experiments led to a “collec¬ 
tive guidance” hypothesis |6], in which a cluster of cells 
where each individual cell has no information about the 
gradient may nevertheless move directionally. In a sense 
that will become clear, cell-cell interactions allow for a 
measurement of the gradient across the entire cluster, as 
opposed to across a single cell. 

In this paper, we develop a quantitative model that 
embodies the collective guidance hypothesis. Our model 
is based on modulation of the well-known contact inhi¬ 
bition of locomotion (CIL) interaction [TUHH], in which 
cells move away from neighboring cells. We propose that 
individual cells measure the local signal concentration 
and adjust their CIL strength accordingly; the cluster 


moves directionally due to the spatial bias in the cell-cell 
interaction. We discuss the suitability of this approach 
for explaining current experiments, and provide exper¬ 
imental criteria to distinguish between chemotaxis via 
collective guidance and other mechanisms where clus¬ 
ters could gain improvement over single-cell migration 
[Hin]. These results may have relevance to collective 
cancer motility m, as recent data suggest that tumor 
cell clusters are particularly effective metastatic agents 
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FIG. 1. Signal-dependent contact inhibition of loco¬ 
motion creates directed motion, a, Schematic picture of 
model and origin of directed motion. Cell polarities are biased 
away from the cluster toward the direction q* = hy 

contact inhibition of locomotion (CIL); the strength of this 
bias is proportional to the local chemoattractant value S'(r), 
leading to cells being more polarized at higher S. See text 
for details, b, One hundred trajectories of a single cell and 
c, cluster of seven cells. Trajectories are six persistence times 
in length (120 min). Scalebar is one cell diameter. The gra¬ 
dient strength is |VS| = 0.025 in these simulations, with the 
gradient in the x direction. 


We consider a cluster of cells exposed to a chemical gra¬ 
dient S'(r). We use a two-dimensional stochastic particle 
model to describe cells, giving each cell i a position r* and 
a polarity pL The cell polarity indicates its direction and 
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propulsion strength: an isolated cell with polarity p* has 
velocity p*. The cell’s motion is overdamped, so the ve¬ 
locity of the cell is p* plus the total physical force other 
cells exert on it, ■ Biochemical interaction be¬ 

tween cells alter a cell’s polarity pb Our model is then: 

dtv^ = p* + ^ (1) 

3¥=i 

= + (2) 

where F'-l are the intercellular forces of cell-cell adhesion 
and volume exclusion, and are Gaussian Langevin 
noises with S(t — f), where the 

Greek indices /x, v run over the dimensions x, y. The 
first two terms on the right of Eq. are a standard 
Ornstein-Uhlenbeck model [HI HH]: p* relaxes to zero 
with a timescale r, but is driven away from zero by the 
noise ^(t). This corresponds with a cell that is orienta- 
tionally persistent over a time of r. 

We have introduced the last term on the right of Eq. 
to describe contact inhibition of locomotion (GIL). GIL 
is a well-known property of many cell types in which cells 
polarize away from cell-cell contact [nilllllMI]- We 
model GIL by biasing p® away from nearby cells, toward 
q® = ^ where f®-^ = (r® — r^)/|r® — r^ | is the unit 

vector pointing from cell j to cell i and the sum over 
j ^ i indicates the sum over the neighbors of i (those 
cells within a distance Dq = 1.2 cell diameters). While 
this is motivated by GIL in neural crest, it is also a natu¬ 
ral minimal model under the assumption that cells know 
nothing about their neighbors other than their direction 
f®-^. For cells along the cluster edge, the GIL bias q® 
points outward from the cluster, but for interior cells q® 
is smaller or zero (Fig. ED- This is consistent with experi¬ 
mental observations that edge cells have a strong outward 
polarity, while interior cells have weaker protrusions [3]. 

Chemotaxis arises in our model if the chemoattrac¬ 
tant S{r) changes a cell’s susceptibility to GIL, /3®, ,0® = 
,5iS'(r®). This models the result of [3] that the chemoat¬ 
tractant Sdfl stabilizes protrusions induced by GIL [3]. 
We also assume that the cell’s chemotactic receptors are 
not close to saturation - i.e. the response is perfectly lin¬ 
ear. If GIL is present even in the absence of chemoattrac¬ 
tant {S = 0), as in neural crest [3], i.e. (3^ = Po + PS{r), 
this will not significantly change our analysis. Similar re¬ 
sults can also be obtained if all protrusions are stabilized 
by Sdfl (t regulated by S), though with some complica¬ 
tions {Appendix, Fig. Al). 

Analytic predictions for cluster velocity. -Out model 
predicts that while single cells do not chemotax, clusters 
as small as two cells will, consistent with [3]. We can 
analytically predict the mean drift of a cluster of cells 
obeying Eqs. m 

(3) 


where the approximation is true for shallow gradients, 
S'(r) « So + r • VS. (• • •)c indicates an average over the 
fluctuating p® but with a fixed configuration of cells r®. 
The matrix Ai only depends on the cells’ configuration, 

^ ( 4 ) 

i 

where, as above, q® = ■ Eq. ^ resembles the 

equation of motion for an arbitrarily shaped object in 
a low Reynolds number fluid under a constant force 
/3rVS [22]: by analogy, we call A4 the “mobility matrix.” 
There is, however, no fluctuation-dissipation relationship 
as there would be in equilibrium [33] . 

To derive Eq. we note that Eq. states that, in 
our units, the velocity of a single cell is equal to the 
force on it (i.e. the mobility is one). For a cluster of N 
cells, the mean velocity of the cluster is 1/iV times the 
total force on the cluster. As F®-^ = —F-^®, the cluster 
velocity is V = ■ When the cluster config¬ 

uration changes slowly over the timescale r, Eq. can 
be treated as an Ornstein-Uhlenbeck equation with an 
effectively time-independent bias from GIL. The mean 
polarity is then (p®) = /3V ■ f®U with Gaussian fluc¬ 

tuations away from the mean, ((p), — (pp)^) = cr^r. The 
mean cell cluster velocity is 

= ( 5 ) 

i j~i 

In a constant chemoattractant field, S = Sq, no net mo¬ 
tion is observed, as linear or slowly- 

varying gradients S{r) k Sq + r ■ VS, and we get Eq. 

Cluster motion and chemotactic efficiency depend on 
cluster size, shape, and orientation.- Within our model, 
a cluster’s motion can be highly anisotropic. Consider a 
pair of cells separated by unit distance along (cos 0, sin 6*). 
Then by Eq. Mxx = \cos^9, Mxy = Myx = 
I cos 9 sin 9, AAyy = | sin^ 9. If the gradient is in the x di¬ 
rection, then {Vx)c = ^ cos^ 9 and {Vy)c = ^ cos0sin6>, 
where Vb = idT|VS'|. Cell pairs move toward higher 
chemoattractant, but their motion is along the pair axis, 
leading to a transient bias in the y direction before the 
cell pair reorients due to fluctuations in p® (Fig. [^. We 
compare our theory for the motility of rigid cell clusters 
(Eq.§ with a simulation of Eq. [T|^ with strongly adher¬ 
ent cell pairs with excellent agreement (Fig.[^. 

For the simulations in Fig.j^and throughout the paper, 
we solve the model equations Eqs. [T]|^ numerically using 
a standard Euler-Maruyama scheme. We choose units 
such that the equilibrium cell-cell separation (roughly 
20 /iin for neural crest 0 ) is unity, and the relaxation 
time T = 1 (we estimate r = 20 minutes in neural crest 
0 )- Within these units, neural crest cell velocities are 
on the order of 1, so we choose <7 = 1- this corresponds 
to a root mean square speed of an isolated cell being 


(V)c « PtM ■ VS' 
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(|V|2)l/2 = 2V2 CTT^/^ ft! 1.4 microns/minute. The typ¬ 
ical cluster velocity scale is Vq = ;0t|VS'|, which is 0.5 
(0.5 microns/minute in physical units) if IVS"! = 0.025 
and S'(O) = 1, corresponding to /3* changing by 2.5% 
across a single cell at the origin. Cell-cell forces are 
chosen to be stiff springs so that clusters are effectively 
rigid (see Appendix for details). 



FIG. 2. Adherent pairs of cells undergo highly 
anisotropic chemotaxis. The average chemotactic veloc¬ 
ity of a highly adherent cell pair (14)c depends strongly on 
the angle 9 between the cell-cell axis and the chemotactic 
gradient. Cell pairs also drift perpendicular to the gradient, 
yf 0. Vb = /3 t|VS| is the velocity scale; IVSI = 0.025. 
Simulations are of Eqs. [l|^ We compute (Vb)c by track¬ 
ing the instantaneous angle, then averaging over all veloci¬ 
ties within the appropriate angle bin. Error bars here and 
throughout are one standard deviation of the mean, calcu¬ 
lated from a bootstrap. Over n = 13, 000 trajectories of Or 
(120 minutes) are simulated. 

We can also compute M. and hence (14) for larger clus¬ 
ters (Table SI, Appendix, Fig. A2). For a cluster of Q 
layers of cells surrounding a center cell, = f{Q)5^ui 
with f{Q) — 2 +^Q+gQ 2 ■ A cluster with Q layers has N = 
1 -I- 3Q -I- 2>Q^ cells; thus the mean velocity of a Q-layer 
cluster is given by (14)/Vb = M = — , where 

M. = \ {Mxx + Aiyy) is the angular average of M. We 
predict that (14)/Vb first increases with N, then slowly 
saturates to 3/2. This is confirmed by simulations of the 
full model (Fig. [^). We note that (14) is an average 
over time, and hence orientation (see below, Appendix). 
We can see why (14) saturates as iV —>• oo by consider¬ 
ing a large circular cluster of radius R. Here, we expect 
q® = an on the outside edge, where a is a geometric 
prefactor and n is the outward normal, with q® = 0 else¬ 
where. Then, ~ (i?d0) n^(6»)ri, = 2ad^j., 

independent of cluster radius R. A related result has 
been found for circular clusters by Malet-Engra et al. 
[3]; we note that they do not consider the behavior of 
single cells or cluster geometry. 

The efficiency of cluster chemotaxis may be measured 
by chemotactic index (Cl), commonly defined as the ratio 
of distance traveled along the gradient (the x displace¬ 
ment) to total distance traveled [24]; Cl ranges from -1 to 
1. We define Cl = (14)/(|V|), where the average is over 
both time and trajectories (and hence over orientation). 
The chemotactic index Cl may also be computed ana¬ 


lytically, and it depends on the variance of V, which is 
{{Vx - {Vx))^) = {(Vy - {Vy))^) = a’^TjN. In our model. 
Cl only depends on the ratio c of mean chemotactic ve¬ 
locity to its standard deviation, 

CI= v^c/Li/2(-cV2) 

( 14 ) PtM\^S\ 

c = , =- . (b) 

where T 1/2 is a generalized Laguerre polynomial. When 
mean cluster velocity is much larger than its fluctuations, 
c 1 and Cl —)■ 1, but when fluctuations are large, 
|c| <C 1 and Cl —)■ 0 {Appendix, Fig. A3). Together, 
Eq.i Eq. [^ and Table SI provide an analytic predic¬ 
tion for cluster velocity and Cl, with excellent agreement 
with simulations (Fig. [^. We note that (14 )/Vq only 
depends on cluster configuration, where Vb = Pt\S/S\, 
so (14(A))/Vb collapses onto a single curve as the gradi¬ 
ent strength is changed (Fig. By contrast, how Cl 

increases with N depends on |V5'| and a (Eq. § Fig.|§3). 



FIG. 3. Larger cell clusters chemotax more effectively, 
but their velocity saturates a, As the number of cells N 
in a cluster increases, the mean velocity {Vx) increases with N 
but then saturates; the mean velocity can be collapsed onto 
a single curve by rescaling by Vb = /3 t|VS'|. b, The chemo¬ 
tactic index GI also saturates to its maximum value. Black 
squares and lines are the orientationally-averaged drift veloc¬ 
ity computed for rigid clusters by Eq. and Eq. [^ Golored 
symbols are full model simulations with strong adhesion. Gell 
cluster shape may influence (W) {Appendix Fig. A4); our cal¬ 
culations are for the shapes in Table SI. Error bars here are 
symbol size or smaller; n > 2000 trajectories of 6r are used 
for each point. 


In our model, clusters can in principle develop a spon¬ 
taneous rotation, but in practice this effect is small, and 
absent for symmetric clusters (see Appendix). 

Motion in non-rigid elusters.- While we studied near- 
rigid clusters above, our results hold qualitatively for 
clusters that are loosely adherent and may rearrange. 
Cell rearrangements are common in many collective cell 
motions pSfBS] , but we note that in |3| clusters are more 
rigid. We choose cell-cell forces F®-! to allow clusters to 
rearrange (see Appendix, |29|), and simulate Eqs. l][2 As 
in rigid clusters, (14) increases and saturates, while Cl 
increases toward unity, though more slowly than a rigid 
cluster (Fig. &b). Clusters may fragment; with increas¬ 
ing X, /3® increases and the cluster breaks up (Fig. [^). 
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Cluster breakup can limit guidance - if /? is too large, 
clusters are not stable, and will not chemotax. 

In Fig. I^b, we compute Cl and velocity by averag¬ 
ing over all cells, not merely those that are connected. 
If we track cells ejected from the cluster, they have an 
apparent Cl > 0, as they are preferentially ejected from 
the high-/3® cluster edge (Appendix). Experimental anal¬ 
ysis of dissociating clusters may therefore not be straight¬ 
forward. Anisotropic chemotaxis is present in non-rigid 
pairs, though lessened because our non-rigid pairs rotate 
quickly with respect to r (Appendix). 





0 200 400 0 200 400 0 200 400 0 200 400 

X [microns] 


FIG. 4. Nonrigid clusters may also chemotax via 
collective guidance, a, As the number of cells A in a cluster 
increases, the mean velocity (14) increases with N but then 
saturates, b, Chemotactic index also approaches unity, but 
slower than in a rigid cluster. Rigid cluster theory assumes 
the same cluster geometries as in Fig. Averages in a-b are 
over n > 20 trajectories (ranging from n = 20 for N = 217 
to n = 4000 for N = 1,2), over the time 12.5r to 50r. c, 
Breakdown of a cluster as it moves up the chemoattractant 
gradient. X marks the initial cluster center of mass, O the 
current center. IVS"! = 0.1, d = 1 in this simulation. 


Distinguishing between potential collective chemotaxis 
models. -Our model explains how chemotaxis can emerge 
from interactions of non-chemotaxing cells. However, 
other possibilities exist for enhancement of chemotaxis 
in clusters. Coburn et al. showed that in contact- 
based models, a few chemotactic cells can direct many 
non-chemotactic ones [14]. If single cells are weakly 
chemotactic, cell-cell interactions could amplify this re¬ 
sponse or average out fluctuations m- How can we 
distinguish these options? In lymphocytes |1], the mo¬ 
tion of single cells oppositely to the cluster immediately 
rules out simple averaging or amplification of single cell 
bias. More generally, the scaling of collective chemo¬ 
taxis with cluster size does not allow easy discrimina¬ 
tion. In Fig. at large N, (14) and Cl saturate. As 
an alternate theory, suppose each cell chemotaxes nois¬ 
ily, e.g. p* = poVS” -I- A*, where A are independent 
zero-mean noises. In this case, (V) = poVS indepen¬ 
dent of N, and ((V^ — (V^))"^) ~ 1/Af, as in our large-A 


asymptotic results and the related circular-cluster the¬ 
ory of |4|. Instead, we propose that orientation effects 
in small clusters are a good test of emergent chemotaxis. 
In particular, studying cell pairs as in Fig. is critical: 
anisotropic chemotaxis is a generic sign of cluster-level 
gradient sensing. Even beyond our model, chemotactic 
drift is anisotropic for almost all mechanisms where single 
cells do not chemotax, because two cells separated per¬ 
pendicular to the gradient sense the same concentration. 
This leads to anisotropic chemotaxis unless cells integrate 
information over times much larger than the pair’s reori¬ 
entation time. By contrast, the simple model with single 
cell chemotaxis above leads to isotropic chemotaxis of 
pairs. 

How well does our model fit current experiments? We 
find increasing cluster size increases cluster velocity and 
chemotactic index. This is consistent with |3], who see 
a large increase in taxis from small clusters (< 20 cells) 
to large, but not |3], who find that Cl is similar between 
small and large clusters, and note no large variations in 
velocity. This suggests that the minimal version of collec¬ 
tive guidance as developed here can create chemotaxis, 
but does not fully explain the experiments of [3- There 
are a number of directions for improvement. More quan¬ 
titative comparisons could be made by detailed measure¬ 
ment of single-cell statistics iniisn], leading to nonlinear 
or anisotropic terms in Eq. Our description of CIL 
has also assumed, for simplicity, that both cell front and 
back are inhibitory; other possibilities may alter collec¬ 
tive cell motion [20]. We could also add adaptation as in 
the LEGI model [STJ |32] to enable clusters to adapt their 
response to a value independent of the mean chemoat¬ 
tractant concentration. We will treat extensions of this 
model elsewhere; our focus here is on the simplest possi¬ 
ble results. 

In summary, we provide a simple, quantitative model 
that embodies a minimal version of the collective guid¬ 
ance hypothesis Eli and provides a plausible initial 
model for collective chemotaxis when single cells do not 
chemotax. Our work allows us to make an unambigu¬ 
ous and testable prediction for emergent collective guid¬ 
ance: pairs of cells will develop anisotropic chemotaxis. 
Although there has been considerable effort devoted to 
models of collective motility [mEUI], ours is the first 
model of how collective chemotaxis can emerge from sin¬ 
gle non-gradient-sensing cells via collective guidance and 
regulation of CIL. 

BAC appreciates helpful discussions with Albert Bae 
and Monica Skoge. This work was supported by NIH 
Grant No. POl GM078586, NSE Grant No. DMS 
1309542, and by the Genter for Theoretical Biologi¬ 
cal Physics. BAC was supported by NIH Grant No. 
P32GM110983. 
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APPENDIX 

CLUSTER CHEMOTAXIS WHEN CHEMOATTRACTANT REGULATES CELL PERSISTENCE r 


Within the main paper, we have assnmed that the chemoattractant concentration S'(r) regulates the susceptibility 
of a cell to contact inhibition of locomotion /3*, with = ,0S'(r*). This models the stabilization of protrusions 
induced by contact interactions. This is consistent with the results of Theveneau et al. [3], who find that protrusion 
stabilization is stronger in clusters than in single cells. However, very similar results can be found if we assume that 
/3 is constant and the signal regulates the time required for the cell’s polarity to relax, i.e. r® = fS{r). In this case, 
the mean polarity of a cell is (p®) = /3t® 

(V) « PfM - VS (t regulation) (Al) 


where the mobility matrix A4 is the same as in the main paper, However, because r varies 

over space, the fluctuations will also vary: {{V^ — {V^))^) = = a'^N~^fS, where S = is 

the mean signal across the cluster. For this reason, the chemotactic index in the r-regulation model will depend on 
VS/S , and will not be constant over a linear gradient. 

In addition, a single cell with a persistence time r that depends on the chemoattractant level will undergo biased 
motion. This is shown in Fig. |Al| below. This drift can be made smaller than the CIL-driven cluster drift, as it is 
independent of /3, while the cluster drift is proportional to /3. 



FIG. Al. Single cells in a spatially-varying r develop a mean drift. The mean x and y velocities for a cell with spatially 
varying r are shown: r = r (So + |VS|a:), with r = 1, So = 1, |VS| = 0.025. Result is average over n = 10® iterations, each 
started at the origin; error bars indicate ([V),(^) “ (W(^))]^)^®^^/\Ai-- 


DERIVATION OF THE MOBILITY MATRICES FOR Q-LAYER OLIGOMERS 


We can compute the mobility matrix of the Q-layer oligomers for arbitrary Q. Our mobility matrix is given by Eq. 
with q® = To simplify the calculation, we can make a few assumptions. First, we note that Aixx = -^yyi 

but Mxy = -Myx = 0 for the Q-layer oligomer. We only need to calculate M = ^ {-Mxx + M.yy). The only cells i 
in the sum of Eq. that are nonzero are those around the boun dary. Ai does not depend on orientation, so we can 
compute the sum oligomer (Eig. A2), then multiply by six. However, this double-counts 

the corner cells, so we must weight them by 1/2. We then hnd M = + §<3] /N{Q), where N{Q) = l + SQ + SQ"^ 

is the number of cells in the cluster. 


We present the mobility matrices for both Q-layer oligomers and other cluster shapes in Table XT 
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FIG. A2. Geometry of Q-layer oligomer, illustrated for Q = 3. The top face is highlighted by a dashed line. 


Shape 

M 

Angularly averaged Al 

Dimer 

oo 

fr:) 

1/4 

Trimer 

OO 

o 

/l/2 0\ 

V 0 1/2^ 

1/2 

Tetramer 

# 

/l/2 0\ 

V 0 3/4^ 

5/8 

Heptamer 

OO 

ooo 

oo 

/6/7 0 \ 

V 0 6/7^1 

6/7 

Q-layer oligomer 

ooo 0=2 

ooo^ 

m° 

/f(Q) 0 \ 
V 0 f(Q)/ 

x(r)\ — 9Q^+3Q 

J — 2+6Q+6Q2 


TABLE Al. Mobility matrices M for several cell configurations. For each of the configurations shown, nearest-neighbor 
cells have unit separation. A Q-layer oligomer has N{Q) = 1-1- 3Q -I- 3Q^ cells. M is given for the orientation shown in the left 
column; other orientations may be found by transforming the mobility tensor; Al = | {Mxx + Myy) (see Section | below). 


ROTATIONAL TRANSFORMATION AND AVERAGING OF THE MOBILITY MATRIX 


We can compute the mobility matrix of a rotated cluster of arbitrary shape from Eq. If we rotate our cluster, 
which we assume is centered at the origin, by an angle 9, (r®) = 7?.(0)r*, we find that 


M' = n^^{e)Mo.pn,p{e) 


(A2) 


where we have assumed the Einstein summation convention and TZ(9) is the rotation matrix ( sin0 j 

y sin0 COS0 j 

matrix terms, Al' = Ti{9) ■ Al • [Tl{d)]^. If we average over 9, we find 


27r 


d9Ai'{9) = - ( -^vv -M-yx 

2 I Myx - M XV M XX + M vv 


(A3) 


We can show from the definition Eq. [^that Al^,, = Al,,^, so the off-diagonal entries of the averaged matrix are zero, 
and therefore A d9M'^j^{9) = \{Mxx + M.yy)5y,v In other words, when averaged over orientation, a cell cluster’s 
mobility matrix is just the constant Al times the identity. 
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COMPUTING THE CHEMOTACTIC INDEX 

We showed in the main paper that within our model, assuming that the cluster rearrangement is slow with respect 
to the polarity dynamics and thus each cell’s polarity is given by a biased Ornstein-Uhlenbeck process, the velocity 
of a rigid cell cluster is 


V = (V) + A 


(A4) 


where A is a Gaussian random variable with zero mean and variance (A^Aj^) = We want to compute the 

chemotactic index. Cl; assuming the gradient is increasing in the x direction, this is 


CI = 


{V.) 

(|V|) 


(AS) 


where the average is both over time and over many trajectories. We note that this is a useful definition for us because, 
in our model, neither (V) nor A depend on the absolute value of the chemoattractant S. More care must be taken 
in other cases. To compute Cl, we need to compute (|V|). |V| is, in our case, given by a Rice distribution, and this 
moment can be calculated. 


(|V|) 



-(A^ + Ag) 
2r2 

{Vy)?] 


(A6) 

(A7) 

(A8) 


where V = + Vy. We now switch to polar coordinates, Vx = V cos(/), Vy = V sin()), and correspondingly write 

(14) = VCOS0 and (Vy) = vsinO, where = (14)^ + Thus, 


(|V|) = 


27rr2 


dcj) / dVV^ exp 
do 


1 


= f2/ dVV^^^T^ 


_^ (1/2 

2r2 


{1/2 + - 21/z/cos(6» - ((i)} 

U.2) Io{Vj^/r^ 


(A9) 

(AlO) 


where Io{x) is the modified Bessel function of the first kind. This integral may be evaluated, resulting in 


(|v|) = ry^Li/2(-i.2/2r2) 


(All) 


where the generalized Laguerre polynomial T 1/2 is given by 

Li/ 2 {x) = e“4 [(1 _ x)Iq{-x/ 2) - xIi{-x/2)] . (A12) 

Within our average over trajectories, we are averaging over the orientation of the cluster; thus we expect {Vy) = 0 
for a chemoattractant gradient in the x direction, and v = (14) = ^rMdxS. T^ = ((I4 — (14))^) = ((14 ~ (14))^) = 
o^t/N. This leads to the result stated in the main paper. 


CI= V2Ac/Ti/2(-c72) 


c = 


(14) 


PrMdxS 


V((^M-(i4)f) 


(A13) 


where in our notation, we could also write c = v/T. We plot the result of Eq. in Fig. A3 below; we see that 
Cl —>■ 1 as c ^ 1 (corresponding to cluster velocities much larger than the noise in cluster velocity) and Cl —>■ 0 if 
|c| <C 1 (cluster velocity much smaller than the noise). We also note that chemotaxis could oppose the direction of 
the gradient (chemorepulsion) - in this case, CI(—c) = —CI(c). 
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FIG. A3. Chemotactic index Cl as a function of the parameter c. 


VELOCITY AND Cl OF IRREGULAR CLUSTERS 


In the main paper, we presented results on the velocity and chemotactic index of Q-layer oligomers. Here, we show 
the velocity and chemotactic index of imperfect clusters. We begin with a Q-lsyev oligomer, and then remove n cells 
at random from the outer layer; this process is repeated 200 times for each n from 1 to 6Q (the number of cells in 
the outer layer). An example is presented in Fig. |A4[ with Q = 5 and n = 5 cells removed. The mobility matrix is 
computed for each cluster, and used to compute (14) and Cl (Fig. A4). We see that though different configurations 
can lead to different mean velocities for the same number of cells, the general trend is captured by the results for 
intact oligomers (dashed line and square symbols in Fig. A4). 
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FIG. A4. Cluster shape effects, in addition to cell number, can affect velocity and Cl. Top: illustration of Q-layer 
oligomer with a few cells removed from the external layer. Bottom: Velocity and chemotactic index for clusters of different 
shapes. Different colors indicate the size of the base cluster from which cells are removed. Black squares connected by dashed 
lines show the results for intact oligomers. For the GI plot, we apply our usual parameters and |VS| = 0.025. All results in 
this figure are theoretical results for rigid clusters only, not full simulations. 
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TRANSIENT ROTATION OF CLUSTERS 

Though we have primarily focused on the translational motion of the cluster, rotational motion can also occur in our 
model, both through rotational diffusion and biased motion. We note that transient rotational events are observed in 
[1]. Under assumptions similar to our main results, clusters have mean angular velocity ~ A • V^, where A depends 
on cluster geometry. This is again similar to an oddly-shaped particle sedimenting in a low Reynolds number flow 
[22] . However, the symmetric clusters in Table SI have A = 0 and do not rotate. If A yf 0, clusters rotate to a fixed 
angle to the gradient direction; there is no persistent rotation in a linear gradient (rotational motion may also be 
suppressed if (3 is large). However, in nonlinear gradients, persistent rotation of asymmetric clusters may be induced. 

We can analyze potential biases by determining the net “torque” Lz = ^ P*] ^ applied to the cluster by the 

cells. This torque is, on average. 


(L,)=^;3r5(U)[dUxq*]^ (A14) 

i 

where q* = and Sr^ = r* — rcm is the displacement from the cluster center of mass. 

What torque is required to cause the cluster to move at a fixed angular velocity? For a cluster moving in a 
rigid rotation with angular velocity H, the cell velocities are v* = H (—Jr*, i5rj,). To achieve this, each cell must 
have a polarity of p* = H (—Jr*, JrJ,), leading to Lz = - |Jr*p. The angular velocity is thus related to Lz by 

LI = Lz/ I We thus find, for linear gradients, ^ -I- r • VS”, 

(H) = ^tA ■ WS (A15) 


where the vector A only depends on the cluster geometry, 

^ ^r* [Jr* X q*] 


(A16) 


where q* is defined as above. (Note that x q* = r* x q* = — 


iriOrj'i = 0, allowing us to drop a center 
A = 0. Cell clusters must lack an inversion symmetry to be 


of mass term.) For all of the shapes listed in Table A1 
rotated by the gradient. 

However, even if A yf 0, clusters will not persistently rotate. We can see that if we rotate the cell cluster around 
its center of mass, A must also rotate as a vector. If th e gra dient is along the x direction, this lets us write 

calculated for a reference geometry. We see that if 


A16 


{LI) = (9) = /3t [^3,(0) cos9 — Ay(0) sind], where A(0) is Eq. 

A yf 0, the cluster will rotate to a stable angle 9* given by tand* = Ax{9i)/Ay{Q). In a linear gradient, there is no 
persistent rotation, though in nonlinear gradients, persistent rotation of asymmetric clusters may be induced. 

We note that Eq. |A15| is not as quantitatively accurate as the corresponding result for translational motion, at least 
for the parameter set in the main paper; this occurs because a small deviation from the equilibrium polarity (p*) can 
create a relatively large change in torque, which often resists rotation. For instance, if the cluster is rotated by a small 


angle J without a corresponding change in (p*), there will be a restoring torque proportional to /3rJ. Thus, Eq. A15 


will be more accurate for systems where the relaxation time r is smaller compared to the rotational timescale of the 
cluster. For similar reasons, for our rigid cluster parameter set, rotational diffusion can be quite slow (Movie SI). 


NONRIGID CLUSTER SIMULATIONS 


In this section, we present additional results on nonrigid clusters. We observe that the cluster anisotropy observed 
in Fig. 2 of the main paper persists in fluid clusters, but is somewhat weaker (Fig. A5). This occurs because the 


rotational diffusion of pairs of cells at our nonrigid parameter set is significantly faster than that of the rigid parameter 
set (compare Movie SI, Movie S2). As the cluster’s polarities (p*) are influenced by the cluster orientation over a 
timescale r, as diffusion becomes faster with respect to r, anisotropy decreases. We also note that we would not 
expect our rigid-cluster results to be precisely accurate even without this effect, as the pair’s separation will fluctuate 
around its equilibrium value, changing Ai. 

We presented in the main paper a figure showing how the cluster velocity and chemotactic index depend on the 
cluster size (Fig. 4 of the main paper). However, in Fig. 4, we treat all of the cells that were initially in contact 
as a single cluster, even if they broke apart - thus we are plotting velocity and Cl versus the total number of cells 
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FIG. A5. Anisotropy exists but is lower in nonrigid clusters. A pair of cells with our nonrigid cluster parameters has 
anisotropic chemotaxis, but with a slightly weaker anisotropy. Trajectories here are measured over the time range 6.25r to 50t; 
n = 20, 200 trajectories are simulated. Rigid cluster results are computed using a cell-cell spacing of 0.75 to roughly match the 
separation seen in simulations, so the curves plotted are 0.75 x | cos^ 6 and 0.75 x |cos0sin0. 


in the simulation. An alternate way to compute a curve showing velocity and Cl as a function of cluster size would 
be to look at a simulation in which an initially large cluster breaks into smaller clusters (as seen in 0), and then 
track the smaller clusters. Doing this can yield different results, because the history of the smaller clusters matters. 
In particular, clusters are more likely to break off from the side of the big cluster at higher S (see e.g. Movie S2) - 
leading to important biases. For instance, we note in Fig. A 61 that smaller clusters (which have been ejected) have a 
larger velocity than large ones, even though isolated small clusters are slower than isolated large clusters. Similarly, 
we see in Fig. |A6[ 3 that even isolated cells develop an apparent chemotactic index. This occurs even though single 
isolated cells in our model have a behavior that is completely independent of the chemotactic signal - if a single cell 
is isolated for a long enough time, its dynamics will again be unbiased. 



FIG. A6. Analyzing clusters in the process of breaking down can lead to apparent single-cell chemotaxis. Tracking 
clusters as they break off from a larger cluster leads to an increased velocity for smaller clusters, and a non-zero chemotactic 
index even for single cells. This figure is generated by simulating clusters initially of 127 cells (6-layer oligomers) over a time of 
50r (1000 minutes), and then computing the average velocity and average instantaneous chemotactic index 14/|V| as a function 
of the number of cells in a given sub-cluster; 200 different trajectories are used. We discard the first 25r of the trajectory, so we 
are focusing on the late cluster breakup stage. Two cells are considered to be in the same cluster if they are within a distance 
Do of one another (e.g. the force between them is nonzero). We only show points in this figure where we have at least 25t of 
trajectory time with clusters of that size totaled over all 200 simulations. Glusters of, e.g. 50 cells are not typically seen, and 
therefore not shown in this figure. 
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NUMERICAL DETAILS OF FULL MODEL SIMULATION 

For our simulations, we solve the model equations numerically using a standard Euler-Maruyama scheme. For our 
rigid cluster simulations, we adapt the cell-cell force from [35] 

f Vr - 1) , < 1 

F*' = , 1 < d*^ < Do (A17) 

[ 0 dd > Do 

where d*-^ = |r* — |. This force is a repulsive spring below the equilibrium separation (which is one in our units, 20 

microns in physical units), an attractive spring above it, and vanishes above Do- Do = 1.2 in all of our simulations, 
and we use Vr = Va = 500. This keeps the clusters very rigid. 

For our non-rigid cluster simulations, we adapt the many-body force chosen by Warren to simulate vapor-liquid 
coexistence [53], 


= [Awid^^/Do) + B{p^ + p>) w{d^^)\ 


(A18) 


where w{r) = (1 — r) for r < 1 and 0 otherwise. The densities are defined by p* = ^jWp{d^^), where the sum is 
over all cells, including i, and Wp{r) = ^(1 — r)^ for r < 1 and 0 otherwise. The force of Eq. A18 is composed of 
an attractive force that goes to zero at a separation of Dq, and a repulsive force that is zero beyond the distance of 
cell-cell overlap (1 in our units). Both attractive and repulsive forces have a finite value even with cells completely 
overlapping (“soft cores”). The strength of the repulsive force increases with increasing cell density - this makes 
the force explicitly dependent on many-body interactions. This force makes developing fluid droplets relatively easy, 
even with short-range interactions [53]. We use A = —23.1 and B = 7.35 for our non-rigid cluster simulations unless 
otherwise noted. 

We initialize our clusters centered at the origin. For rigid clusters, we start our simulations with the shapes given 


in Table A1 but rotated to a random angle, and a spacing of the equilibrium spacing (unity) for rigid clusters. For 


non-rigid clusters, we start with the appropriate Q-layer cluster at a random angle, but with a spacing of 0.7. For 
non-rigid clusters, we initialize 2-, 3-, and 4-cell clusters by removing the appropriate number of cells randomly from 
the outer layer of a heptamer. In both cases, we initialize the polarity p to a random value from its distribution (i.e. 
/3Vq® plus an appropriate noise). 


TABLE OF PARAMETERS 


Parameter symbol 

Name 

Value in our units 

T 

Persistence time 

1 

O 

Characteristic cell speed (OU noise parameter) 

1 

P 

CIL strength 

20 for rigid clusters, 1 for nonrigid clusters 

Va 

Adhesion strength 

500 (rigid clusters only) 

Vr 

Cell repulsion strength 

500 (rigid clusters only) 

Do 

Maximum interaction length 

1.2 

So 

Signal strength at origin 

1 

At 

Time step 

10““^ for rigid clusters, 5 x 10“^ for nonrigid clusters 

A 

Attraction strength for Warren force (Eq. jAlSj) 

-23.1 (nonrigid clusters only) 

B 

Repulsion strength for Warren force (Eq. |A18|) 

7.35 (nonrigid clusters only) 


TABLE A2. Parameters used 
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